
#ifndef _UPWIND_XLD_H_
#define _UPWIND_XLD_H_

#include "Equations/Momentum_XCoDi.hpp"

namespace Tuna {
  
  template <typename T, int Dim> class Upwind_XLD;
  
  template <typename T>
  struct Typeinfo <Upwind_XLD<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };
  
  template <typename T>
  struct Typeinfo <Upwind_XLD<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };
  
  template <typename T>
  struct Typeinfo <Upwind_XLD<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };

template<class Tprec, int Dim>
class Upwind_XLD : public Momentum_XCoDi<Upwind_XLD<Tprec, Dim> > 
{
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::aE;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::aW;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::aN;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::aS;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::aF;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::aB;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::aP;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::sp;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::dx;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::dy;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::dz;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::bi;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::ei;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::bj;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::ej;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::bk;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::ek;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::dt;    
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::applyBoundaryConditions1D;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::applyBoundaryConditions2D;
    using GeneralEquation<Momentum_XCoDi<Upwind_XLD<Tprec, Dim > > >::applyBoundaryConditions3D;
	
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::Gamma;
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::Rayleigh;
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::phi_0;
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::T;
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::u;
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::v;
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::w;
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::p;
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::du; 
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::calc_du_2D;
    using Momentum_XCoDi<Upwind_XLD<Tprec, Dim > >::calc_du_3D;

public:

  typedef Tprec prec_t;
  typedef typename TunaArray<prec_t, Dim >::huge ScalarField;

  Upwind_XLD() : Momentum_XCoDi<Upwind_XLD<prec_t, Dim> >() { }
  ~Upwind_XLD() { }; 

  inline bool calcCoefficients1D() { };
  inline bool calcCoefficients2D();
  inline bool calcCoefficients3D() { };
  inline void printInfo() { std::cout << " Upwind_XLD "; }
};


template<class Tprec, int Dim>
inline bool Upwind_XLD<Tprec, Dim>::calcCoefficients2D()
{
    prec_t dy_dx = Gamma * dy / dx;
    prec_t dx_dy = Gamma * dx / dy;
    prec_t dxy_dt = dx * dy / dt;
    prec_t ce, cw;
    prec_t cn, cs;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	{
	    ce = ( u(i+1, j) + u(i,j) ) * 0.5 * dy;
	    cw = ( u(i-1, j) + u(i,j) ) * 0.5 * dy;
	    cn = ( v(i,j) + v(i+1,j) ) * 0.5 * dx;
	    cs = ( v(i,j-1) + v(i+1,j-1) ) * 0.5 * dx;

	    if ( ce > 0 ) ce = 0.0; 
	    else          ce = -ce;
//
// This statement:	    
	    if ( cw <= 0 ) cw = 0.0;
//
// is more efficient than the next similar one:
//	    if ( cw > 0 ) cw = cw; 
//	    else          cw = 0.0;

	    if ( cn > 0 ) cn = 0.0;
	    else          cn = -cn;
	    if ( cs <= 0 ) cs = 0.0; 

	    aE (i,j) = (dy_dx + ce);
	    aW (i,j) = (dy_dx + cw);
	    aN (i,j) = (dx_dy + cn);
	    aS (i,j) = (dx_dy + cs);
	    aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS (i,j) + dxy_dt;
//		+ (ce - cw) + (cn - cs);
// Term (ce - cw) is part of discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.
	    sp (i,j) = u(i,j) * dxy_dt - ( p(i+1,j) - p(i,j) ) * dy;
	}    
    calc_du_2D();
    applyBoundaryConditions2D();
    return 0;
}

//
//---------------------  3D  ---------------------
//
/* *
template<class Tprec, int Dim>
inline bool Upwind_XLD<Tprec, Dim>::calcCoefficients3D() 
{
    prec_t dyz = dy * dz, dyz_dx = Gamma * dyz / dx;
    prec_t dxz = dx * dz, dxz_dy = Gamma * dxz / dy;
    prec_t dxy = dx * dy, dxy_dz = Gamma * dxy / dz;
    prec_t dxyz_dt = dx * dy * dz / dt;
    prec_t ce, cw;
    prec_t cn, cs;
    prec_t cf, cb;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int k = bk; k <= ek; ++k)
	for (int i =  bi; i <= ei; ++i)
	    for (int j = bj; j <= ej; ++j)
	    {
		ce = ( u(i+1, j, k) + u(i,j,k) ) * 0.5 * dyz;
		cw = ( u(i-1, j, k) + u(i,j,k) ) * 0.5 * dyz;
		cn = ( v(i,j,k) + v(i+1,j,k) ) * 0.5 * dxz;
		cs = ( v(i,j-1,k) + v(i+1,j-1,k) ) * 0.5 * dxz;
		cf = ( w(i,j,k) + w(i+1,j,k) ) * 0.5 * dxy;
		cb = ( w(i,j,k-1) + w(i+1,j,k-1) ) * 0.5 * dxy;
				
		if ( ce > 0 ) ce = 0.0; 
		else          ce = -ce;
//
// This statement:	    
		if ( cw <= 0 ) cw = 0.0; 
//
// is more efficient than the next similar one:
//	    if ( cw > 0 ) cw = cw; 
//	    else          cw = 0.0;
		if ( cn > 0 ) cn = 0.0;
		else          cn = -cn;
		if ( cs <= 0 ) cs = 0.0; 
		if ( cf > 0 ) cf = 0.0;
		else          cf = -cf;
		if ( cb <= 0 ) cb = 0.0;

		aE (i,j,k) = (dyz_dx + ce);
		aW (i,j,k) = (dyz_dx + cw);
		aN (i,j,k) = (dxz_dy + cn);
		aS (i,j,k) = (dxz_dy + cs);
		aF (i,j,k) = (dxy_dz + cf);
		aB (i,j,k) = (dxy_dz + cb);
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) + aN (i,j,k) + aS (i,j,k)
		    + aF (i,j,k) + aB (i,j,k) + dxyz_dt;
//		+ (ce - cw) + (cn - cs) + (cn - cs);
// Term (ce - cw) is part of discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.
		sp(i,j,k) += u(i,j,k) * dxyz_dt - 
		  ( p(i+1,j,k) - p(i,j,k) ) * dyz; 
	}    
    calc_du_3D();
    applyBoundaryConditions3D();   
    return 0;
}
/* */

} // Tuna namespace

#endif //_UPWIND_XLD_H_

















